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TEMP02, a new pulsar timing package. Ill: Gravitational 
wave simulation 
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ABSTRACT 

Analysis of pulsar timing data-sets may provide the first direct detection of gravi- 
tational waves. This paper, the third in a series describing the mathematical framework 
implemented into the tempo2 pulsar timing package, reports on using tempo2 to sim- 
ulate the timing residuals induced by gravitational waves. The tempo2 simulations 
can be used to provide upper bounds on the amplitude of an isotropic, stochastic, grav- 
itational wave background in our Galaxy and to determine the sensitivity of a given 
pulsar timing experiment to individual, supermassive, binary black hole systems. 

Key words: methods: numerical - gravitational waves - pulsars: general 



1 INTRODUCTION 

Sazhin (1978) and Detweiler (1979) were the first to realise 
that pulsar timing observations provide a powerful tool for 
detecting ultra-low frequency (f g ~ 10 -9 Hz) gravitational 
waves (GWs). The precision with which millisecond pulsars 
are now being timed makes it possible that pulsar timing 
experiments could provide the first direct detection of a GW 
signal 1 . The Parkes Pulsar Timing Array (PPTA) project 
(e.g. Hobbs 2008, Manchester 2008 and references therein) is 
an attempt to achieve this ambitious goal by making regular 
observations of 20 bright millisecond pulsars. 

Recent theoretical work (e.g. Jaffe & Backer 2003, 
Wyithe & Loeb 2003) suggests that the strongest signal 
potentially detectable by such experiments would be an 
isotropic stochastic GW background caused by coalescing 
supermassive black holes in the centres of merging galaxies. 
Jenet et al. (2005) showed that in order to detect this signal, 
the 20 PPTA pulsars will need to be timed to a precision of 
~0.1 /xs over a timespan of ~5 yr. To date, the PPTA project 



1 Observations of the first binary pulsar, B1913+16 (Hulse & 
Taylor 1974), provided the first evidence for the existence of GW 
emission. The pulsar timing experiments described in this paper 
are designed to make a direct detection of GWs. 



has data spanning ~3 yr with root-mean-square (rms) resid- 
uals of typically 0.1— 3/is, but it is expected that these resid- 
uals will significantly improve over the next few years with 
new observing systems and enhanced signal processing pro- 
cedures. Therefore, it is now appropriate to determine how 
these existing data-sets can be used to limit the amplitude 
of GW signals and how future data-sets will be analysed in 
order to detect a GW signal and determine its properties. 

Pulsar observations lead to measurements of pulse 
times-of-arrival (TOAs; ^ bs ) at an observatory. Paper I 
(Hobbs, Edwards & Manchester 2006) and Paper II (Ed- 
wards, Hobbs & Manchester 2006) of this series detail how 
the new pulsar timing package, tempo2^ 
ta bs to the proper time of emission, t^ ST , 



, is used to convert 
as 



(1) 



A© is the transformation required to convert the site ar- 
rival times to the solar system barycentre, Ars is the excess 
propagation delay due to the interstellar medium and Ab is 
the transformation to the pulsar frame for binary pulsars. 
Tempo2 compares the derived time of emission with a pulsar 



2 The TEMP02 software and documentation are available from our 
website http: //www. atnf . csiro . au/research/pulsar/tempo2. 
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model to form "timing residuals" , which are equivalently the 
deviations between the observed TOAs and the model pre- 
dictions. For a perfect pulsar model, random receiver noise 
and no other systematic effects, these timing residuals will 
have a mean of zero and be uncorrelated, corresponding to 
a flat, or "white", spectrum. Since tempo2 does not include 
GW sources in the timing model, the existence of any such 
sources will induce a signal in the timing residuals. The aim 
of this paper is to describe how this signal can be simulated 
and how such simulations aid searches for GW signals within 
our existing data-sets. 

Since the intrinsic pulsar pulse period, spin-down, or- 
bital motion and various astrometric parameters are a priori 
unknown, they must be determined from the pulsar tim- 
ing data. In common with other pulsar timing analysis pro- 
grams, tempo2 uses initial estimates of the pulsar param- 
eters to obtain "pre-fit" timing residuals and then uses a 
least-squares fitting procedure to fit an analytical model to 
obtain improved pulsar parameter estimations and "post- 
fit" timing residuals (full details are given in Paper I). The 
net outcome of this process is that a polynomial and vari- 
ous spectral components are removed from the post-fit tim- 
ing residuals. Any GW signal with a period larger than the 
time-span of the data is largely absorbed by the removal of 
the low-order polynomial terms. Hence, pulsar timing ex- 
periments are only sensitive to GW signals with periods less 
than, or equal to, the time-span of the data (typically years), 
corresponding to frequencies in the range 1-30 nHz. 

The three basic types of GW sources that have been dis- 
cussed in the literature are (1) continuous wave sources (Pe- 
ters 1964), (2) burst sources (e.g. Thorne & Braginskii 1976, 
Damour & Vilenkin 2001, Kocsis et al. 2006 and Enoki & 
Nagashima 2007) and (3) stochastic backgrounds (e.g. Jaffc 
& Backer 2003, Wyithe & Loeb 2003, Maggiore 2000). The 
GW strain spectrum for a stochastic background is thought 
to be a power- law in the GW frequency, f g , as 

hc{h) = Aa {jt) ' (2) 

where /i yr = 1/lyr and A g is dimensionless. For a back- 
ground generated by supermassive binary black holes, a = 
-2/3 and A g ~ 10~ 15 (Jaffe & Backer 2003, Wyithe & Loeb 
2003). Standard models of inflation (e.g., Turner 1997; Boyle 
& Buonanno 2007) produce GW backgrounds with ampli- 
tudes well below detectable limits with current experiments 
(A g ~ 10~ 18 ), but some non-standard models (e.g., Gr- 
ishchuk 2005) have a ~ — 1 and A g ~ 10~ 15 . Cosmic string 
cusps are also expected to produce a GW background with 
a = —7/6 and A g can become as large as 10 -14 (Damour & 
Vilenkin 2001; Caldwell, Battye & Shellard 1996). 

Determining a rigorous limit on A g is not trivial as 
real pulsar data-sets have irregular sampling, non-white 
noise due to instrumental problems, intrinsic pulsar timing 
noise, astrometric and orbital parameter fitting, and inac- 
curacies in the terrestrial time-standard or in the planetary 
ephemeris. Jenet et al. (2006) recently described how simu- 
lating GW signals within tempo2 allows rigorous limits to 
be placed on A g , which take into account the majority of the 
issues affecting real pulsar observations. The use of tempo2 
and the methods employed were only outlined in the Jenet el 
al. paper; full details are provided here. Unfortunately, the 
Jenet et al. technique can only be applied to timing residu- 



als that have a white spectrum. We have recently developed 
a new technique that makes no assumption on the spectrum 
of the timing residuals. This recent work will be presented 
in a forthcoming paper. 

In §2 we provide the mathematical framework that 
allows tempo2 to simulate GW sources. This is divided 
into sections considering the timing effects induced by non- 
evolving GW sources (§2.1) and evolving sources (§2.2). In 
§3 we demonstrate applications of this mathematical frame- 
work within tempo2. 



2 SIMULATING THE EFFECT OF GW 
SOURCES ON PULSAR TOAS 

The equations presented in this paper describe how the in- 
duced timing residuals for a given pulsar due to a GW signal 
can be calculated. However, this is not sufficient for our pur- 
poses. We must be able to simulate the effects of a GW on 
the actual pulse TOAs, because the process of fitting a tim- 
ing model to obtain the residuals will modify the effects of 
a GW. Numerous methods exist within tempo2 to simulate 
such TOAs. These methods are all based on the following it- 
erative procedure. First, a set of observation dates and times 
are defined by the user. Second, the entire tempo2 timing 
procedure (as described in Paper II) is carried out in or- 
der to obtain pre-fit timing residuals. This procedure uses a 
user-specified timing model defining the pulsar being simu- 
lated and assumes that the dates and times described above 
represent pulse TOAs. Third, these pre-fit timing residuals, 
which really describe the timing model, are subtracted from 
the original arrival times. The goal is to obtain arrival times 
which, when fitted with a timing model, give zero residu- 
als. However, because of various non-linear operations in the 
modelling and fitting process, this procedure must be iter- 
ated until the resulting pre-fit residuals are adequately close 
to zero for the simulation being planned. These TOAs can 
subsequently be modified by the addition of white Gaussian 
noise, a model of the pulsar timing noise and/or the GW 
signal. The final TOAs are stored as if they were actual pul- 
sar observations and can be processed using standard fitting 
and analysis routines. 

tempo2 employs two different techniques to simulate 
the effect of GWs on pulsar timing residuals. The first tech- 
nique is used for constant frequency (i.e. "non-evolving") 
sources. The second is used for simulating the effects of bi- 
nary systems that are evolving. Since the latter technique is 
computationally expensive, the former is used to calculate 
the effects of a stochastic background of GWs. The basic 
implementation into tempo2 is described below. We also 
provide detailed derivations of all the main equations in the 
Appendix. 

2.1 Non-evolving GW sources 

The majority of the GW sources that may be detectable by 
pulsar timing are expected to evolve over timescales much 
longer than the typical observation time. Hence, the non- 
evolving algorithm used in tempo2 can be used in most 
GW simulations. 

The non-evolving GW simulation algorithm in tempo2 
has been defined so that the user can input pulsar and GW 
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Figure 1. Configuration of the coordinate system used through- 
out this paper. Note that a X S = f. 



source positions in an equatorial coordinate system. Figure 1 
represents the position, with respect to the Earth, of either 
a pulsar (with unit position vector f p and distance d v ) or a 
GW source (f 9 , d g ). The sources are specified by their right 
ascension and declination (a;, Si). 

Tempo2 assumes a globally flat coordinate system with 
three spatial coordinates (x,y,z) and one temporal coordi- 
nate t. GWs are treated as a tensor field in this background 
space-time. A single plane GW takes the form 



him = Re 



(3) 



where the indices I and m range from 1-3 corresponding 
to the three spatial coordinates. Ai m is a constant tensor 
amplitude, k g is the three-dimensional GW vector and co g is 
the GW angular frequency. A GW signal causes fluctuations 
in the observed pulsar's observed spin frequency Sf/f. The 
induced pulsar timing residuals are given by the integral of 
this quantity over time. The timing residuals induced by a 
GW of the above form are given by 



R(t) 



-±Re 



c 



, (4) 



where Q = 1 — cos 6 and 6 is the angle between the direc- 
tion of the pulsar and the direction of the GW source (see 
Appendix) . 

In tempo2, the GW tensor amplitude is specified in the 



coordinate system where the GW is propagat- 



ing along the —r g direction. GWs consistent with Einstein's 
equations have two independent degrees of freedom, which 
we label as A+ and A x . Written in terms of these values, 
the tensor amplitude takes the form: 



Since tempo2 allows one to arbitrarily specify the entire ten- 
sor amplitude, one can generate GWs consistent with any 
general metric theory. Once the GW amplitude is specified in 



the (r g ,&g,5 g ) coordinate system 



is evaluated by 



transforming both f p and Ai m into the global (x,y,z) coor- 
dinate system. This scalar quantity is then used in equation 
4 to calculate the induced pulsar timing residuals for the 
given pulsar. In the remainder of this section, we will dis- 
cuss how the above general framework is used to simulate 
GWs from a single, non-evolving, binary system as well as 
from a stochastic background of GWs. 

2.1.1 GWs from supermassive black-hole binary systems 

Supermassive black-hole binary systems in the cores of 
galaxies are expected to be sources of detectable GWs. For 
long-period binary systems, the time it takes for the orbital 
period to evolve under the action of GW emission (~ 10 4 
years for a system with chirp mass M c — 10 9 solar masses 3 
and a three-year orbital period) is much longer than any 
reasonable observation time. Hence, the binary system may 
be treated as non-evolving. In general, the GWs emitted by 
a binary system will be elliptically polarised (Blanchet et al. 
1996). Since the tensor amplitude is a complex quantity, the 
effects of such GWs can be calculated using the framework 
described above. 

In the current tempo2 implementation, only systems 
with zero eccentricity are considered for the non-evolving 
case. This is a valid assumption since binary systems tend to- 
wards zero eccentricity much faster then the decay timescale 
(Peters 1964). 

Following Wahlquist (1987), TEMP02 models GWs 
emitted from a binary system by setting A + and A x as 
follows: 



A x = 
where 



,-ie n ^ 3 + cog cos ( 2 0-, + i4 cos (0.) S in(20)] (6) 
e -»9„ + CQS 0.) sin ( 2 0) _ i4 COS (6i) cos(20)] (7) 



Ag = 



(8) 



Oi is the orbital inclination angle, <j> is the orientation of the 
line of nodes, 8 n is the orbital phase angle at the line of 
nodes, M c is the binary chirp mass, lj is the orbital fre- 
quency and d g is the distance to the source. Note that the 
GW angular frequency, uj 3 = 2lo . 

2.1.2 A stochastic background of GWs 

It is possible within tempo2 to specify a large number of 
individual GW sources, each with different properties. A 
stochastic background of GWs is simulated by randomly 
specifying the source directions and tensor amplitudes of 
the GWs generated by these sources. Such a background is 
described by its characteristic strain spectrum, h c (f) (equa- 
tion 2). In order to simulate such a background, probability 
distributions for the GW parameters are defined as follows. 



\ / \ 3 / 5 

A+ A x (5) 3 The chirp mass is defined as M c = (mi +m 2 ) ^ J^^p j 

A x — A+ / where mi and m,2 are the masses of the binary companions. 
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The source directions are chosen uniformly on the celestial 
sphere so that the respective probability distribution func- 
tions are given by: 

1 

2' 



P(sin5) 



(9) 
(10) 



The GW frequencies are chosen to be uniformly distributed 
in logtJ 9 : 



"a log(^) 





UJl < LUg < UJh 

otherwise 



(11) 



where uii and u h can be defined by the user, but default to 
tu h — 27r/(ld) and u>; = 27r2^ where T is the time-span of 
the observations. 

A+ and A x , the parameters used to determine the ten- 
sor amplitude (see equation 5), are treated as real numbers 
(i.e. the imaginary parts are set to zero) that are normally 
distributed with zero mean and rms given by 



\og(u> h /u)i) 



N 



M/) 



(12) 



where N is the number of individual plane waves used to 
generate the background. 

Given the above choice of distributions, the simulated 
background will be isotropic, unpolarised and have a Gaus- 
sian amplitude distribution with characteristic strain h c (f). 
For the existing simulations within TEMP02, h c (f) is taken 
to be of the form given by equation 2. The spectral index, 
a, and the amplitude, A g , depend on the physical processes 
generating the background and may be specified by the user. 



2.2 Evolving sources 

In general, a binary system will evolve under the action 
of GW emission and can have non-zero eccentricity. Al- 
though it is possible to use the above framework to model 
such a source, it is not very convenient. A separate mod- 
ule (the GWevolve plug-in; see below) has been developed 
in tempo2 to deal with this case. Full details of the equa- 
tions integrated numerically within TEMP02 were provided 
by Jenet et al. (2004) and therefore are not reproduced here. 
As shown in §3.2 the user inputs the initial eccentricity and 
orbital periods to obtain the resulting timing residuals gen- 
erated using the specified geometry of the orbit. 

Unfortunately, solving the differential equations nu- 
merically is computationally expensive. Hence, this module 
is currently only used to simulate well-defined, individual 
evolving sources. 



3 APPLICATIONS WITHIN TEMP02 

As described in Paper I, the tempo2 software is based 
around 'plug-ins' that add to the functionality of the pack- 
age. The mathematical framework described above allows 
for the development of new plug-ins to simulate, study, and 
detect GW signals. A listing of the current plug-ins available 
for GW research is presented in §3.4. These existing plug- 
ins are divided into those 1) simulating the induced timing 




Frequency (day 1 ) 

Figure 2. Average power spectrum obtained from 1000 re- 
alisations of white timing residuals and fitted using the 
PSR J0437— 4715 timing model. The vertical dotted lines corre- 
spond to periodicities of lyr and 0.5 yr respectively. The dashed 
lines correspond to the orbital period of 5.7 d and twice the orbital 
period respectively. 



residuals due to single GW sources or from a stochastic GW 
background, 2) producing an upper bound on the amplitude 
of any GW background, 3) determining the sensitivity of a 
given set of pulsar timing residuals to single GW sources 
and 4) for inspecting the resulting timing residuals. In this 
section, we demonstrate these plug-ins. 

It is important that tempo2 is used for fitting the pul- 
sar's astrometric, pulse and, if applicable, orbital parameters 
when studying the induced timing residuals due to GW sig- 
nals as such fits reduce the detection sensitivity at various 
characteristic frequencies. In Figure 2 we show the average 
power spectrum obtained after fitting a standard pulsar tim- 
ing model to a white-noise, daily sampled data-set with an 
rms timing residual of 100 ns. The absorption features are 
due to the removal of power by fitting for the astrometric, 
rotational and orbital parameters. It is difficult to obtain 
a straight-forward description of these spectral features as 
they depend on the details of the fitting procedure and on 
the sampling of the data. However, in the tempo2 routines 
for simulating and studying GW signals that are described 
below, detailed analytic descriptions of such spectral fea- 
tures are not required; the simulations take all such effects 
into account. 



3.1 Stochastic GW background 

GW backgrounds can be simulated using the GWbkgrd 
plug- in. The power-law spectrum of h c (f) leads to a power- 
law spectrum for the pulsar timing residuals with spec- 
tral exponent a ros = 2a — 3 (see equation A58). Hence, 
for a background generated by supermassive black holes 
where a = —2/3, the induced timing residuals will have 
a much steeper red-noise spectrum with spectral expo- 
nent a ICS = —13/3. Example timing residuals (simulated 
every two weeks for 3000 d) are shown in Figure 3 for 
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J1939+2134 (rms = 0.614 /xs) post-fit 




-1000 1000 

MJD-51501.2 



Figure 3. Example timing residuals induced by a stochastic 
GW background (with A g = 10 — 14 ) after fitting for the pulsar's 
(PSR B1937+21) pulse period and its first derivative. The error 
bars correspond to 100 ns of additional white, Gaussian noise. 

PSR B1937+21 after fitting for the pulsar's pulse frequency 
and its first derivative. 

In Figure 4 we show the power spectrum 4 of 512 weekly 
sampled simulated residuals induced by a GW stochastic 
background with A g — 10 -15 and a — —2/3. The simulation 
was repeated 1000 times and the average power spectrum is 
shown, with the theoretical spectrum (^4g/127r 2 )/~ 13/ ' 3 yr 3 
drawn as a solid line. In each simulation 10000 plane GWs 
were summed as discussed in §2.1.2. On this scale the av- 
erage power spectrum can barely be distinguished from the 
theoretical line, except at high frequencies. The apparent 
high frequency noise corresponds to an rms of 0.2 ns and 
occurs due to rounding errors in the pulsar timing model 
computations. Since tempo2 has been designed to maintain 
1 ns precision, and our best observations currently have an 
rms residual of 50 ns, this white noise is negligible. 

3.1.1 Producing an upper-bound on the background 

Many techniques have been described in the literature for 
determining an upper bound on A g . The earliest work 
(e.g. Kaspi et al. 1994, McHugh et al. 1996) was based on 
analysing the measured post-fit timing residuals. More re- 
cently, Jenet et al. (2006) used the tempo2 simulations of 
GW backgrounds that are described in this paper to pro- 
duce an upper bound on A g for various values of a. This 
technique has limitations. Notably it requires that the ob- 
served timing residuals are 'white' (defined as being a data- 

4 Analysis of such steep 'red' spectra is challenging because of the 
irregular sampling of the observations and spectral leakage from 
the low-frequency components. In this case the sampling is reg- 
ular and leakage was eliminated by prewhitening the time series 
with a second difference filter and postdarkening the spectrum 
with the inverse of the transfer function of the second difference 
filter. We have found that most observations can be handled with 
combinations of interpolation and prewhitening. These techniques 
are being integrated into TEMP02 and will be discussed elsewhere. 
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Figure 4. The average spectrum of 1000 GW background real- 
isations for A g = 10 — 15 and a = —2/3 for 512 weekly-spaced 
simulated observations. The solid, diagonal line is the theoretical 
spectral density. 

set whose power spectrum is independent of frequency, or 
equivalently, for which the data points have no temporal 
correlation). Each data-set used in the Jenet et al. (2006) 
work was tested by 1) forming power spectra (constructed 
using a Lomb-Scargle periodogram and using Gram-Schmidt 
orthonormal polynomials) and searching for significant pe- 
riodicities and 2) averaging adjacent points to confirm that 
the variance of the timing residual decreases with the num- 
ber of points averaged. However, the power spectrum at low 
frequencies is suppressed by the fitting procedure carried out 
by tempo2 and therefore even though a data-set may pass 
the tests described above, it may not have a purely white 
spectrum. 

For completeness, we describe here the details of the 
tempo2 usage in the original Jenet et al. (2006) method, 
but emphasise that new techniques are currently being de- 
veloped that are not restricted to white data-sets. It is ex- 
pected that an implementation of many of these new tech- 
niques (e.g. van Haasteren et al., in press, Anholm et al., in 
press) will use and develop the tempo2 functionality that is 
described below. 

In the Jenet et al. (2006) method, a statistic is first 
defined that is sensitive to a GW background. Following 
the terminology of the original paper we define each pulsar 
data-set to consist of n p measured residuals, x p (i), a time 
tag t p (i) and an uncertainty o~ P (i) where i is the data sample 
index and p is an index referring to a particular pulsar. Each 
data-set may be unevenly sampled. Normalised time tags 

r p (i) = 2(t p (i) - Cn)/(C ax - CiJ - 1 (13) 

are defined where t?. is the earliest time and t^ax the time 
of the most recent observation for pulsar p. Hence, r p (i) 
runs from —1 to 1. These r p (i) values are used in a weighted 
Gram-Schmidt orthogonalisation procedure to determine a 
set of orthonormal polynomials, j l p (i), defined from 

n \ 

#(»)#(<) , n4 s 

i=Q 

where j l p (i) is the Z'th order polynomial evaluated at r p (i) 
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and Sik is the standard Kronecker delta function. The follow- 
ing coefficients are calculated using the orthonormal poly- 
nomials, j l p{i), and the timing residuals, x p (i): 

n p — l , 

C » - 2^ alii) ■ (15) 

i— 

The pulsar average polynomial spectrum is given by 

P^y-Kt (16) 

V 

where the weighted variance, v p , is defined as 
A- X^T-o 1 ( a; p(*) — ^'p) 2 /°pM an d ^ i s ^ ne mean °f x. For a 
stochastic background dominated by low-frequency noise, Pi 
will be large for low values of I. Hence, T = X);=o ^ * s use d 
as a statistic to detect the background. The upper limit, n, 
can be selected by the user, but n = 7 was used throughout 
the Jenet et al. (2006) paper. 

The background will be "detected" if T > T where 
To is set so that the false-alarm probability is given by Vf. 
By default Vf =0.1%. To is obtained using the following 
Monte-Carlo procedure. First, standard pulsar timing proce- 
dures are followed to obtain "pre-fit" timing residuals, R\{i), 
for each pulsar data-set. These are subtracted from the origi- 
nal site-arrival-times and the procedure iterated until arrival 
times, t\(i), are obtained that are exactly predicted by each 
pulsar's timing model. Noise is then added back to the ar- 
rival times. Since only pulsar residuals that are consistent 
with "white noise" can be analysed by the Jenet et al. (2006) 
method, an independent data set with the same noise dis- 
tribution as the original is obtained by adding a shuffled 
version of the timing residuals R\(i) to t\{i). With this new 
simulated set of site-arrival-times, the entire tempo2 tim- 
ing procedure is repeated in order to obtain a new set of 
"post-fit" timing residuals, i?f(i). The detection algorithm 
is subsequently applied to -Rf(i) and the output statistic 
Tj is recorded. This procedure is repeated for Nit iterations 
where Na is set, by default, to 10000. These Tj values are 
subsequently inspected to determine To. 

Finally, the upper bound on A g is determined so that 
the probability of detecting the background with A g — 
Supper is Vd- By default, Vd = 95%. This upper bound is 
determined by adding a GW background of a given ampli- 
tude to if(i). As above, the fitting procedures are carried 
out to obtain "post-fit" timing residuals and the detection 
algorithm applied to obtain T. If T > To then the back- 
ground has been detected. The amplitude is changed using 
a bracketing procedure in order to determine the amplitude 
Supper which gives a detection probability of Va- 
in this technique, the actual timing residuals are used 
only as a mechanism for generating instances of white noise 
in the simulations. If the spectrum of the measured timing 
residuals is red then this technique will provide an upper 
bound which is too low because the shuffled observations 
(which will be white) will give lower detection statistics than 
a simulation based on the correct noise spectrum. A plug-in 
package, CHECkWhite, is available in tempo2 to test the 
"whiteness" of a data-set; see §3.4. 

The default values of Na and N gv , have been chosen to 
produce a stable upper limit that has the precision neces- 
sary for current astrophysical applications. To demonstrate 



this, we use the data set for PSR J1857+0943 that was first 
described by Kaspi et al. (1994) and used to determine an 
upper bound by Jenet et al. (2006) of A < 1.45 x 10~ 14 (cor- 
responding to a bound on the energy density per unit loga- 
rithmic frequency interval of f2 gw [l/(8yr)]/i 2 < 1.3 x 10~ 7 ) 
for a = — 1. Multiple simulations using the same observa- 
tions, but with different realisations of the GW background 
and with different shuffles of the data, produces a mean up- 
per bound of A < 1.54 x 10~ 14 and standard deviation of 
0.06 x 10 -14 . It should be noted that the "whiteness" of 
the residuals of PSR J1857+0943 is suspect because the ob- 
served detection statistic is 2.4 times higher than the mean 
of the simulated detection statistics using shuffled observa- 
tions. A detection statistic would exceed this value only 3% 
of time by chance, suggesting that the residuals are some- 
what red. 

3.1.2 Detecting the background 

Hellings & Downs (1983) showed that a GW background sig- 
nal can be detected by searching for correlations in the tim- 
ing residuals of many pulsars spread over the sky. Within the 
framework of general relativity, the induced timing residuals 
for any isotropic, stochastic GW background are correlated 
with a well-defined zero-lag angular correlation function: 

c{fi) = \x\ixx-\ + \ + \s[x) (17) 

where x = [1 — cos#]/2 for angle 9 on the sky between two 
pulsars. Our simulations successfully reproduce this angu- 
lar correlation. In Figure 5 we show the results of simu- 
lated timing residuals (using the GWbkgrd plug-in) in the 
presence of a GW background for the 20 PPTA millisecond 
pulsars (no pulsar noise is added). For each pulsar pair, we 
plot the zero-lag correlation versus the angular separation 
of the pulsars on the sky. The solid line is the predicted 
functional form (equation 17). In order to produce the fig- 
ure using standard correlation techniques, we have selected 
the GW spectral exponent a = +3/2 which corresponds to 
arcs = (i.e. white timing residuals) and have simulated 
regularly sampled timing residuals with two weekly sam- 
pling over five years. In general, obtaining the zero-lag cor- 
relations between pulsar pairs in the presence of red-noise 
in the timing residuals and uneven sampling is challenging 
and requires pre-whitening of the data in order to attain 
the maximum possible signal-to-noise ratio. Hence, Figure 5 
displays the optimal effect that could be achieved with pre- 
whitening techniques. These algorithms, and their implica- 
tions for GW background detection, will be described in a 
subsequent paper. 

3.2 Simulating single sources and the effect of 
parameter fitting 

Tempo2 plug-ins are available to simulate both non-evolving 
individual GW sources (GWsingle) and evolving sources 
(GWevolve). The non-evolving source simulations can eas- 
ily be shown to produce sinusoidal residuals of the correct 
amplitude and phase for a given GW source and pulsar posi- 
tion. In Figure 6 we use the GWevolve plug-in to reproduce 
the expected PSR B1855+09 timing residuals for the pos- 
tulated binary supermassive black-hole system in the radio 
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Figure 5. Pair-wise angular correlation curves for 20 simulated 
pulsar data-sets in the presence of a gravitational wave back- 
ground with power-law index a = +3/2 and amplitude A g = 0.01. 



galaxy 3C66B (Sudou et al. 2003). As described in Jenet et 
al. (2004) , the induced signal has a low-frequency component 
due to the GW signal at the pulsar and a higher frequency 
component due to the GW signal at Earth. Figure 6a shows 
the pre-fit timing residuals induced by the simulated GWs. 
Figure 6b gives a realistic representation of observed post- 
fit timing residuals if 3C66B did contain a binary black-hole 
system with the chirp mass and period given in Sudou et 
al. (2003). As concluded by Jenet et al. (2004), such a sig- 
nal would be easily detectable, but has not been observed 
in actual pulsar data-sets (Figure 6c). Note that the low- 
frequency term would be indistinguishable from the cubic 
variations often observed and attributed to pulsar period 
irregularities. 

One of the best known candidates for a supermassive 
binary black hole system emitting GWs with frequencies 
detectable by pulsar timing is in the blazar OJ287, where 
a periodicity of ~ 12 yr has been identified in optical out- 
bursts (e.g., Sillanpaa et al. 1996). The parameters of this 
system are not well-defined. However, to obtain an order- 
of-magnitude estimate of the induced timing residuals due 
to the GW emission from this system, we can use the pa- 
rameters originally suggested by Sillanpaa et al. (1988) and 
make no cosmological corrections. They model the system 
with mi = 2 x 10 7 M Q , m 2 = 5 x 10 9 M , an orbital period of 
9 yr in the rest frame of the blazar and an initial eccentricity 
of e = 0.7. The source has a red-shift of 0.306 corresponding 
to a distance of ~1250Mpc. The GWevolve plug-in shows 
that the induced timing residuals due to this system are 
significantly less than 1 ns and therefore undetectable in all 
existing data-sets. Note that cosmological effects will only 
change this result by about 20%. 



3.3 Public data sets 

Many publications which have described limits on the ex- 
istence of a GW background, or on the existence of indi- 
vidual GW sources, have relied on publically available pul- 
sar timing residuals (in particular, most have used the data 
sets made available by Kaspi et al. 1994). It is likely that 
the resurgence of interest in pulsar timing arrays and GW 



detection will lead to many more techniques being devel- 
oped. In order to aid comparison between different tech- 
niques we have made available a set of simulated pulsar 
timing residuals with and without the addition of a GW 
background. These timing residuals, arrival time files and 
parameter models are available from our web-site . The tim- 
ing residuals for the simulated PPTA data are also available 
as an electronic supplement to this paper (see Appendix B). 

• Simulated PPTA data: based on the design specifica- 
tions for the PPTA project, we provide data-sets with two- 
weekly sampling of 20 pulsars with white, Gaussian noise 
giving 100 ns rms timing residuals. We include data-sets 1) 
without the addition of a GW background, 2) with a back- 
ground where A g = 10~ 14 and a = —2/3 and 3) with a 
background where A g = 10~ 5 . The pulsar parameter files 
were obtained from the ATNF pulsar catalogue (Manchester 
et al., 2005). 

• Simulated global timing array data: we provide data- 
sets which are likely to be created by the global pulsar tim- 
ing array (i.e. combining observations of both Northern and 
Southern hemisphere pulsars). We simulate 30 pulsars (the 
20 PPTA pulsars and 10 more Northern millisecond pul- 
sars), with data spans ranging from 5yr to 12 yr, realistic 
observation dates and rms timing residuals from 100 ns to 
1 fis. We provide two different GW background amplitudes 
(with A g = 10~ 15 and A g = 10~ 14 respectively). 

• Simulated SKA data: The SKA is likely to be able to 
time at least 100 millisecond pulsars with rms timing resid- 
uals around 50 ns. In order to simulate possible data-sets 
we select the 100 fastest recycled pulsars in the ATNF pul- 
sar catalogue that are not associated with globular clusters. 
We simulate weekly sampled, white timing residuals over a 
data-span of 10 yr. 



3.4 Plug-in packages 

The following plug-ins are available for tempo2 from our 
website. 

• fake: As described in §2, this plug-in allows the user 
to simulate pulse arrival-times at an observatory that are 
in accordance with a specified pulsar timing model to bet- 
ter than 1 ns. This plug-in has been used in producing the 
publically available files that are described in §3.3. 

• GWbkgrd: This plug- in allows the user to simulate the 
pre- and post-fit timing residuals resulting from a specified 
GW background. 

• GWsingle: Allows the user to simulate the pre- and 
post-fit timing residuals resulting from a non-evolving super- 
massive black-hole binary system at a given distance. 

• GWevolve: This plug-in determines pulse arrival- 
times that have been affected by a binary source evolving 
due to emission of gravitational radiation. 

• GWwhiteLimit: This plug-in implements the tech- 
nique first used by Jenet et al. (2006) to place an upper- 
bound on the amplitude of a GW background. 



5 select the "publically available data files" link from 
http://www.atnf.csiro.au/research/pulsar/tempo2 
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Figure 6. (a) Simulation of the timing residuals induced in the Kaspi et al. (1994) timing residuals for PSR B1855+09 due to the 
postulated supcrmassivc binary black hole system in 3C66B. (b) Simulated residuals after including the GW signal, the measured timing 
residuals and their uncertainties and fitting for the pulsar's spin-down, astrometric and orbital parameters, (c) The observed timing 
residuals. 



• checkWhite: A plug-in to test the "whiteness" of 
a particular data-set. This plug-in plots various power- 
spectral estimates (including a Lomb-Scargle periodogram 
and a Gram-Schmidt orthogonal polynomial power spec- 
trum) and calculates the statistic used in the Jenet et al. 
(2006) upper-bound technique for the actual timing residu- 
als and for shuffled realisations of the timing residuals. 

• plk: This plug-in is available with the default tempo2 
distribution. It allows the user to view pre- and post-fit pul- 
sar timing residuals. The user may turn on (or off) fitting 
for various model parameters and re-calculate post-fit tim- 
ing residuals. Figures 3 and 6 in this paper were obtained 
using this plug-in. 



4 CONCLUSIONS 

A major advantage of Tempo2 over previous pulsar tim- 
ing packages is that its functionality can be expanded using 
plug-in packages. Numerous plug-ins have now been devel- 
oped in order to simulate and analyse the effects of GW 
signals on pulsar timing data. This code has already been 
used to place the most stringent constraints to date on the 
existence of a GW background (Jenet et al. 2006). It is now 
being used to study how a GW background could be de- 
tected, to determine the sensitivity of a given pulsar timing 
array to single and burst GW sources and to study the pos- 
sibilities of pulsar timing array projects with future instru- 
ments such as the Square Kilometre Array telescope, with 
which we hope to not only detect gravitational waves, but 
also to study their properties in detail and the sources from 
which they emanate. 
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APPENDIX A: MATHEMATICAL DESCRIPTION OF GW SOURCES 

This Appendix shows how the expected power spectrum from a stochastic background of GWs is calculated. This power 
spectrum is related to the characteristic strain spectrum, h c , and the normalized power per logarithmic frequency interval, 
0, gw . In order to establish a consistent, well-defined notation, the calculations are presented from first principles. Note that 
we are using standard geometrised units where c = I . 



Al The stochastic background and its energy density 

GWs are linear perturbations to a background space-time metric. For the purpose of this paper, we will assume that the 
background space-time is flat. Hence, the space-time metric may be written as: 

g^ v = r\^ v + (AI) 
where 

-10 
,0100. 

'""'= I o o i o ! 

1 

and hpt, is a small perturbation. The linearised Einstein equations with r)^ as the background space-time take the following 
form: 

h X \^ v — h^^^v — h\,\/j, + h^v^x — 0. (A3) 

A stochastic background of GWs is made up of a sum of plane waves travelling in several directions. Hence, one can write 
the metric perturbation due to a stochastic background as: 



hf,, v = Re 



3=0 



(A4) 



where N is the total number of GWs and A^ v , , kj and ujj are the complex amplitude, spatial wave vector and angular 
frequency of the jth GW respectively, ujj is taken to be positive. 

The stress-energy tensor for a metric perturbation is given by a 4-D volume average: 

T% = Zs / K„, a h%d 3 xdt, (A5) 

where L and T are the spatial and temporal limits of integration respectively. L and T are taken to be several times the 
longest wavelength involved. Using the metric of a stochastic background, equation A4, the energy density, p gw , takes the 
form: 

P 3W = T 7 = 3—^3 / I E {-i^A^--^ + i^A^e x ( A6) 

ji 

I — loJiA 1 ^ e ' 1 +iljiA^ e 1 ' Id xat. 

As long as there are a finite number of plane GWs, or sources, in the sum, one can make the following approximation with 
reasonable accuracy: 

I_L J e ii ^- R ' ) - s - ( ^-"' )t dtd 3 x = 6 jl (A7) 
where 5ji = 1 if j = I and zero otherwise. Using this, the energy density becomes: 

^ = ^E^-^r- (as) 

i 

A 

Next, the above sum will be written in integral form using a probability density function. The amplitude of a given GW 
depends on k and a set of other parameters denoted as a. Examples of these other parameters are mass and distance. Letting 
dP/d n ad 3 k be the probability density for all the parameters on which a general GW may depend, the ensemble-averaged 
energy density is given by 

p3W = lkl ^Uk,3)A^(k,3)N^-d-ad 3 k. (A9) 

Since d 3 k = uj 2 cLodfl where dQ — dcos6d(f> and 6 and <f> are the usual spherical coordinate angles specifying the GW 
propagation direction, the energy density per unit frequency is given by 



10 
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d P g 



du) 



= _L/ 

64tt J 



oj i A^(k,a)A^(k,a)N 



dP 



d n ad 3 k 



d n adfl. 



(A10) 



A2 The stochastic background and the induced timing residuals 

The action of a GW slightly alters the arrival times of radio pulses emitted by a radio pulsar. Equivalently, the rate of arrival of 
the pulses will fluctuate. Since the action of gravity does not depend on the frequency of the electromagnetic (EM) radiation, 
the problem of determining the change in the rate of arrival of pulses of EM radiation simplifies to the problem of finding the 
change in frequency of a single-frequency plane wave or photon. 

The four-dimensional path of a photon is the shortest path between specified end points. The four-dimensional path is 
represented by x^(X) where A is the so-called "affine parameter". At any given point along the path, the wave four- vector is 
determined by fc£ = dx^/dX. Since the infinitesimal distance between two points on the four-dimensional curve is given by 
y/kpkpg^dX, the "distance" between two points along any such curve is given by: 



An 



dX. 



(All) 



Defining L — y/—k^kpg^ v , the shortest path between two fixed endpoints is given by the four Lagrange equations: 

ddL_dI L=Q (A12) 
dX dk a dx a K ' 

Using 



dh 

dk a 
dh 
dx a 



L 
1 1 

12 



k p g a v 



the Lagrange equations yield: 
dkp a _ 1^ , p ,„ 



(A13) 
(A14) 

(A15) 



In order to calculate the terms due to the action of the GW alone, we will assume that both the pulsar and the observer 
are at rest in the global background coordinate system. In this case, k p o is the frequency of the photon. k p o at the pulsar 
will be written as Lu e while at the receiver it will be denoted as u) r . Next, we will write the metric using equation Al and let 
kp — kp + Sk^ where ftp 1 is the photon four- vector in the unperturbed space-time and Skp 1 is the induced perturbation to the 
path. The equation for the perturbed photon frequency is then given by: 



k k h/j,i/fi. 



dSkpo 
dX 

Using equation A4 for the metric, results in 



1t„t 

2 



d8u> 
~dX 



-Re 



. 3 



-itok l t ( -.k 7i aAu, v je ^ 



where £ M (A) is the unperturbed photon path given by 
x M (A) = kp(X — X e ) + x%, 



(A16) 



(A17) 



(A18) 



x'i is the location of the photon emitter (i.e. the pulsar) and A e is the affine parameter of the emitter. Putting this into 
equation A17 and integrating yields: 



5u) r — 3uj e 



-Re 
2 



£ 



where x% is the location of the receiver. Using 

% r X e — kp [\r " -^e) 

together with the unperturbed light travel time between the pulsar and the receiver, D 



(A19) 



(A20) 

c2, it can be shown that 



(A21) 



With the above, equation A19 may be written as 
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5u) r — 8co e = 



Next, take the receiver location to be xp = (t, 0, 0, 0) and write k^jkp 1 as 

k^jkp = — LO e LOj(l — COsd), 



(A22) 



(A23) 



where 8 is the angle between the direction of the GW source and the pulsar. With the above, the fractional frequency shift 
may be written as 



Suj r — Sui e 1 

= -Re 

0J e 2 



Ei.fi," / 1 „iw,U(l-cose,') 

2 A„^e 



1 — cos t 



(A24) 



The coordinate system has been chosen so that each A^ v only has spatial components. Using this, the final form of the 
fractional frequency shift is: 



8uj r — Sul e 1 

= -Re 

0J e 2 



^ ^ kpkp Ai m jG 



I _ gi^j-Dfl-cosflj) 

1 — cos 9j 



(A25) 



where k p is the unit vector in the direction of the pulsar. 

The change in the arrival time of a pulse at time t is given by the integral of the fractional change in frequency of the 
pulse rate. Hence, the timing residuals induced by a set of plane GWs is given by 



R(t) 



f 

Jo 



Su)r{t') — Side ,., 1 _ 

— at = — Re 

UJe 2 



E 



h l h m A, ■ 



1) 



jjD(l-cos 8j) 



1 — cos ( 



In order to simplify notation, the following definitions are made: 



B{t)j = 



Cj = 



1 i (e 



1) 



a — itjj D(l — cos 0j ) 

1 — cos Oj 



(A26) 



(A27) 

(A28) 
(A29) 



Using the above notation and explicitly taking the real part of the summand, the induced timing residuals become: 
R W = ~\ Yl B ^) C ^ + B*{t)C*E*. (A30) 



Next, the ensemble-averaged power spectrum of the residuals is calculated. For a given length of data, T, the variance of 
the residuals is given by 



± jf R\t)dt- (jf J R(t)d^j 



(A31) 



In order to take the ensemble average of the above, it is assumed that no two GWs have the same W 1 and that GWs with 
different fc M are not related to each other; GWs from different regions of the sky are uncorrelated. Mathematically, the above 
statements are expressed as: 



(Almj Apqk) — 

(AlmjApqk) = (Al m jA pq j)Sjk- 

With the above, the ensemble average variance may be calculated by putting equation A30 into A31: 

2 N 



- 2 ) = k ( \Ci I 2 \Ej | 2 U £ \Bj (t)fdt - i /; Bj {t) 



where \x\ is the complex amplitude of x. The integrals in the summand take the form: 



L \ B] (t)\ 2 dt- ^[ B 3 {t) 



4 % 2 



I — sine' 



Of) 



where sinc(:r) = sin(x)/a;. Using the same technique to derive equation A9, equation A34 may be written as: 



1 — sine' 



(¥) 



\k'k m Ai m (k,3)\ 2 



cos[a;D(l — cos 8)] 
(I - cose) 2 



N 



dP 



d n ad 3 k 



d n ad kdQ,. 



(A32) 
(A33) 



(A34) 



(A35) 



(A36) 
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The fact that k 2 — lo 2 implies d 3 k = uj 2 du>. Hence, the above equation tells us that the power spectrum of the residuals 
is given by 

da 2 



1 



\k l k m Ai m {k, a)\ 2 



1 - cos[wL>(l - cos6l)] dP 



(1 - cose) 2 



d n ctd 3 k 



d n adQ,. 



(A37) 



A3 An isotropic, unpolarised, GW background 



Until now, the derived expressions for both the energy density of a GW background and the induced pulsar timing residuals 
have allowed for an arbitrary directional dependence. Here, the calculations will be simplified for the case of an isotropic 
background and the power spectrum of the induced timing residuals will be calculated in terms of the normalised energy 
density per unit logarithmic frequency interval, Q gw (f). In this case, dP/d n ad 3 k does not depend on the direction of the GW. 
The energy density per unit frequency may be written as: 



cLo 



A* llv A IJ,v can be expressed in terms of the amplitudes of the two independent GW modes A + and A x : 
A;„A»" =2\A + \ 2 +2\A x \ 2 . 

Since the GW background is taken to be unpolarised, (|A + | 2 ) = (\A X | 2 ). Hence 

(a;„a»") = 4\a + \ 2 

and the energy density per unit frequency may be written as 



dpg 



du> 



\A+(u;,a)\ 2 N 



dP 



d n ad 3 k 



d n a. 



(A38) 



(A39) 



(A40) 



(A41) 



In order to calculate the induced residual power spectrum for an isotropic background, equation A37 will be expressed 
in a standard spherical coordinate system with the pulsar located along the z-axis. As before, 9 represents the angle between 
the pulsar and the direction of the GW as well as the standard spherical coordinate polar angle, r is the unit vector pointing 
in the direction of the source. 9 and <j> are unit vectors pointing in the direction of increasing 9 and 4>, respectively. These unit 
vectors, which depend on 9 and (f>, make up a local right-handed coordinate system with 9 x <j> = f. Each Aij can be written 
in terms of the f , 9, <f> coordinate system: 



Agg = —A^ = A4 

= A 60 = A x 



with all other components equal to zero. Using the above, one finds that 

\k l k m A lm (k,d)\ 2 =sin 4 (e)jA+| 2 . 

Since the pulsar lies in the z direction, <j> • z = 0, and 9 ■ z = — sine. The induced timing residuals therefore become: 



da^ 

dio 



2 fuT\ 



\A+\ 2 N 



dP 



d n ad 3 k 



cos[ujD(1 — cose)] 
(1 - cose) 2 



dQ. 



(A42) 
(A43) 
(A44) 

(A45) 
(A46) 



Since the background is assumed to be isotropic, neither dP/d n ad z k nor |A+| 2 depend on direction, hence they are taken 
outside the dfl integral. The integration over solid angle is given by: 



siti4g l-cos[^(l-cose)] ^ 



(1 - cose) 2 



16tt 
3 



8tt 47rsin(2o;D) 



(coD) 2 



(luD) 3 



(A47) 



At this point, the short wavelength approximation will be made (i.e. ujD » 1) so that the last two terms in the above 
are negligible. The power spectrum of the induced timing residuals can now be written as: 

da 2 
dio 



4tv 



\A+\ 2 N 



dP 



d n ad 3 k 



771 

d a. 



(A48) 



Comparing this to equation A41, the power spectrum of the induced residuals may be written in terms of the energy 
density per unit frequency: 

'2 



da' 
dio 

da 2 



I6n 1 dp gw 
3 lo 4 dui 



In terms of frequency (/ = lo/2tt), this becomes: 



1 1 dp g 



df Stv 3 f* df 



smc 



(A49) 



(A50) 
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Typically, the energy density spectrum is written per unit logarithmic frequency interval as: 



da _ 1 1 dp gw 
In terms of 9U) (/) 



1 . 2 ( 2nfT \ 
l_ slnc 



— f, Pa 7f\ i one has 

p c dlog(/) ' 



eicr 



1 — sine' 



(A51) 



(A52) 



where p c — SHq/S-k, and Ho is the Hubble constant. 



A4 Timing residuals and the characteristic strain spectrum 



Several investigators use the "one-sided" strain spectrum, Sh{f), of the GW background or the characteristic strain spectrum 
h c (f). These quantities are defined as: 



f 

Jo 



hc(f) = yffSrfJ). 
Using the same techniques employed above to calculate dp gm /df, one finds that 

= \J A »Aa, k)A^(a, k)N ,„ dP „ , u; 2 divdQd n a 



(A53) 
(A54) 



64tt 4 / f\A+\ 2 N 



dP 



d n ad a k 



d n ad a k 



d n adf 



(A55) 



where the last equality holds for the case of an isotropic, unpolarised background. Using the above and the definition of Sh(f), 
one finds that 



S h lf)=32n 4 f / \A+\ 2 N 



dP 



-d a. 



d n ad 3 k 

Using this together with equation A48, the power spectrum of the residuals is given by: 
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For the case of a power-law characteristic strain spectrum as given by equation 2, the power spectrum of the residuals may 
be written as: 
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Also note that the normalised power per logarithmic frequency interval, Q gw (f), can also be written in terms of Sh(P 
and the characteristic strain spectrum (see equation A52): 
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APPENDIX B: SIMULATED PPTA DATA SETS 

An electronic supplement to this paper includes the simulated PPTA data sets. Three tables are provided giving two-weekly 
sampling of the 20 PPTA pulsars with the addition of 100 ns white, Gaussian noise. The first table has no additional GW 
signal, the second table includes a GW background where A g = 10~ 14 and a = —2/3 and the third table has a background 
where A g — I0~ 15 . The first column in the online tables gives the MJD of the simulated observation and the remaining 20 
columns give the timing residuals for each of the 20 PPTA pulsars in the following order: PSRs J0437— 47f5, J0613— 0200, 
J07I1-6830, J1022+1001, J1024-07f9, J1045-4509, JI600-3053, JI603-7202, J1643-1224, J1713+0747, J1730-2304, 
J1732-5049, J1744-1134, J1824-2452, JI857+0943, J1909-3744, JI939+2I34, J2I24-3358, J2129-5721 and J2145-0750. 



